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Abstract 

The dynamics of two dimensional (2D) vortex fluctuations are investigated 
through simulations of the 2D Coulomb gas model in which vortices are rep- 
resented by soft disks with logarithmic interactions. The simulations strongly 
support a recent suggestion that 2D vortex fluctuations obey an intrinsic 
anomalous dynamics manifested in a long range 1/i-tail in the vortex correla- 
tions. A new non-linear TV-exponent a, which is different from the commonly 
used AHNS exponent clahns and is given by a = 2aAHNS — 3, is confirmed by 
the simulations. The results are discussed in the context of earlier simulations, 

experiments and a phenomenological description. 
PACS numbers: 05.40.+j, 74.40.+k, 74.76.-w, 75.40.Mg 
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I. INTRODUCTION 



Superconducting films, 2D Josephson junctions, and 4 He films undergo a Kosterlitz- 
Thouless(KT) type transition from the superfluid to the normal state. (MB This transition 
is driven by thermally created vortex-antivortex pairs which start to unbind at the transition. 
The high-T c materials can to some extent be regarded as weakly coupled superconducting 
planes, which raises the question to what extent the 2D vortex fluctuations are important 
also for this new class of materials. || In order to assess such questions it is important to 
understand the properties of the vortex fluctuations in the pure 2D case. One motivation 
for the present investigation is to gain such understanding. 

The KT-transition is driven by the vortex fluctuations and this means that the large 
distance and long time behaviour in the transition region is dominated by the properties of 
the vortices. The static properties at the transition are described by the Kosterlitz renor- 
malization group (RG) equations (4| and are rather well understood. |5|,|2| The properties 
related to the dynamics of the vortices constitute a much more open question. We have in 
the present paper performed extensive simulations on a simple dynamical model of vortex 
fluctuations in order to gain some further insight. 

So far the most widespread view on the dynamics of vortex fluctuations derives from series 
of papers by Ambegaokar et al. || We refer to this as the AHNS phenomenology. Somewhat 
later a variant was devised by one of us which we will refer to as the MP-description. || 
These two phenomenologies in fact give very different predictions. Experimental data seem 
to favor the MP-description. [@,[[]|§ Particularly clear experimental evidence for this was 
given by Theron et al in case of a 2D array of Josephson junctions. |§ The MP-description 
has also been clearly borne out in computer simulations on 2D XY-type models. |[T0|JTT[| 
It has been argued that at the heart of the MP-description is a long range 1/t-tail in the 
vortex correlations below the transition temperature. |12[ This 1/t-tail reflects an anomalous 
diffusion of the vortex fluctuations and is then the key to the difference between the AHNS 
and the MP. In the present paper we demonstrate that this anomalous diffusion is already a 
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property of the "simplest" possible dynamical model of vortex fluctuations: A charge neutral 
system of positive and negative particles with logarithmic particle interaction and Langevin 
dynamics. 

Another main prediction of the AHNS is the exponent ciahns of the non-linear IV- 
characteristics i.e. V oc I aAHN $ , However, given the anomalous diffusion and the 1/i-tail 
of the vortex correlations, it has been argued that this prediction is no longer correct. [12 



A scaling argument suggests that the non-linear JV^-exponent a consistent with the 1/t-tail 



is given by a = 2clahns — 3. JEU This prediction is borne out to an excellent degree in the 



present simulations. |13 | 



The content of the present paper is as follows: In section 2 we describe the model and 
the simulation procedure. Sections 3-6 contain the results from the simulations. In section 
3 we show that the model and the simulations correctly reproduce the static KT-transition 
per se. In particular we verify the power law decay of the correlations for large distances 
in the low temperature phase. Good agreement with the predicted power law indices are 
found. Then in section 4 we present the results for the non-linear JV^-exponent a and verify 
the correctness of the new scaling exponent. In section 5 we verify the 1/t tail in the vortex 
correlations. In particular we show how the decay of the temporal correlations depend on 
the wavevector k for small k. Section 6 gives the frequency dependence of the basic linear 
response function describing the coupling to an external electromagnetic field. The results 
are shown to be very well represented by the functional form of the MP description in the 
small k limit. We also show how the response function crosses over to a a more Drude 
like behaviour as the wavevector k is increased. Finally section 7 contains some concluding 
remarks. 

II. MODEL AND SIMULATIONS 

In accordance with AHNS we assume that the dynamics of the vortices to good approx- 
imation is described by the Langevin equation 
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where r is the position of a vortex and F to t is the total force acting on it due to all the 
other particles as well as any externally imposed force, D is the diffusion constant, T is the 
temperature (unit system such that the Boltzmann constant ks = 1), and i] is a random 
force obeying 

(t{t)r] p {t')) = 2DT8 a p8{t - t') (2) 

where a and (3 denote the Cartesian components. This equation describes the strong friction 
limit in which the vortex motion is perpendicular to the applied external current and should 
be a good approximation of a 2D superconductor. |||],[14|,|2| 

According to the vortex-Coulomb gas analogy, the vortices can be described as a gas of 
2D Coulomb charges with logarithmic interaction. |2j The two possible vorticities s ± 1 of a 
2D vortex corresponds to positive and negative Coulomb gas charge. The thermally created 
vortex configurations have zero total vorticity which corresponds to a neutral Coulomb gas. 
H In our model the Coulomb gas charges are taken to be disks of extension r . These disks 
correspond to the vortex cores and are such that the force acting between two particles i 
and j with charges Sj and Sj respectively (in units such that the charge is s = ±1) and 
separated by the distance r is given by 

= s iSj {- - -K x {r/r )) (3) 
r r 

where K\ is a modified Bessel function of order 1. This means that the charge distribution of 
a Coulomb gas particle is soft, which is in accordance with the precise vortex- Coulomb gas 
particle analogy. |l5j Consequently the force between two particles vanishes for r = and is 



proportional to 1/r for r >> r$. Alternatively one may express the two particle interaction 
corresponding to Eq. (|3|) in terms of a potential U (r) 

U{r) = - \n{r/r ) + K {r/r ) (4) 

where K is a modified Bessel function of order 0, 



Fi . = -s.s.blJLufo.) (5) 

' ij ij 

and Yij is the position vector from particle % to particle j. In the present paper length is in 
units of r and time in units of t = Tq/D. 

The simulations are performed for a fixed number of particles N and constant temper- 
ature T. The particles are contained in a 2D quadratic box of side length L with periodic 
boundary conditions. The numerical solutions were obtained by discretizing time into small 
time steps At and introducing a random noise r](t) which acts independently on each particle 
at each time step. The Langevin equation (|1|) is then turned into a finite difference equation 
for the particle system 

N 

Ti(t + At) = n(t) + At ]T F(r y (t)) + AtF ext (t) + ^(t), (6) 

where the indices % and j numerate the particles and the diffusion constant D has been 
absorbed into the time scale and the random force. F(rij(t)) = SiSj—Fij is the force acting 
at time t on particle i due to particle j and F ext is any external force. The random force 
in Eq.(^) can thus be treated as a random displacement vector rji(t) which obeys (compare 
Eq-(D) 

(tf(i)^) = 2r^-0. ( 7 ) 

and is sampled from a Gaussian distribution. This equation is then solved on the computer 
by using a standard Euler integration method. [TO For each temperature of interest the value 
of At was halved repeatedly until no dependence of the time step could be monitored (usually 
At ~ 0.01 to)- The number of time steps needed for convergence is usually 1 — 5 x 10 6 , 
but for the long time correlation functions as many as 15 x 10 6 steps were needed in order 
to obtain decently converged time tails. In practice one has to strike a balance between 
choosing At small enough to ensure that the equation of motion is correctly solved yet as 
large as possible in order to achieve as large time sequences as possible. In practice we 
have been able to meet these conditions below T c without too much problem. However, 



just at and slightly above T c this turned out to be very computer time consuming. Another 
practical problem is to keep track on the influence of the boundary Here particular care has 
to be taken because the two particle interaction is long range i.e. U(r) oc lnr for large r. 
To this end we found it expedient to modify the interaction by a large distance exponential 
cut off A c which could then be varied in order to check the dependence on the largest length 
scales. Thus Eq.(|]) was modified into 

Fij = SiSji^-Kiir/K) - -K^r/ro)) (8) 

corresponding to 

U(r) = K (r/X c ) - ^ (r /A c ) - K (r/r ) (9) 

Typical parameters in the simulations are N = 512 and L/r m 320 which correspond to a 
particle density n = 5 x 10~ 3 r-Q 2 . The ratio X c /L = 0.35 turned out to be an efficient choice. 
The size dependence of the results was checked by varying L for fixed n and ratio X c /L. The 
size L/tq = 320 was in practice large enough to avoid finite size effects except very close 
to the phase transition. In fact we found that simulations onaiV = 512 system were for 
practical purposes large enough for the parameter range we are investigating. However, to 
be on the safe side a fair amount of the numerical data was obtained for iV = 1024 and 
occasional checks for N = 2048 were also performed. 

The Coulomb gas is often discussed in terms of a fugacity variable z where z 2 /A 2 is 
the probability of creating a dipole pair with particle separation ro and A is the phase 
space division for a particle. This means that in our model there exists a non trivial relation 
between nr 2 and z(nr^, T). However, in our present simulations n and T are the fundamental 
variables. 

The basic correlation function which we obtain from the simulations is the Fourier trans- 
form of the charge density correlation function g(r,t) defined as 

1 N 

^tj^E^^ 1 - ( 10 ) 
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In principle this function has a slight directional dependence on k due to our choice of 
periodic boundary on a quadratic box. However, in practice our simulation results are to 
good approximation spherical symmetric so that g(k, t) = g(k,t). 

The results are conveniently discussed in terms of the complex frequency dependent 



dielectric constant 



of the Coulomb gas model which is the basic response function. 



This is related to the correlation function g(k,t) by 
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(13) 



The first equation (|TTD gives the static result which contains the information on the thermo- 
dynamic KT-transition. The following two, ( |l"2"D and (pT3|), contain the information specific 
to the dynamics of the model. 

In addition to the linear response given by Eqs ( ]TT| - |T3"|) we calculate the non-linear re- 
sponse for the case when the model is subject to an external force F ext = SjE where E is 
constant in space and time. In this case we calculate the average particle charge current I p 
per particle 



N 



N [ ^ Sl dt 



1 D 



N 



NT; 



(14) 



i=l 



where the first equality is the definition of I p and the second follows directly from Eq.(p]). 
F^ t is the total force acting on the particle i and the brackets (} denote a time average. 



The results from these simulations are presented in the following three sections. 



III. KT-TRANSITION. 



We will first focus on the static dielectric function l/e(k) = l/i(k,uj = 0) given by 
Eq . (|i~l~l) . This function is related to the linearly screened two particle interaction by 
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1 _ Ueffjk) 

m " ~ow { } 

The "bare" interaction U(k) is in our case given by (compare Eq.flj^)) 

provided L = oo. In practice we use the numerical transform for finite L and periodic 
boundary conditions. The linearly screened interaction is for small k given by [TJJ 

^/W = i FTA^Wj (17) 

where A < A c is the screening length. Consequently we expect that the static dielectric 
function for small k is of the form 

_j_ = i _m_ 

e(Jfe) e /c 2 + A- 2 1 ' 

Figure 1 shows data for l/e(fc) obtained from our simulations for a sequence of temperatures 
at a fixed density n. The filled circles represent the data and the full curves are fits to Eq. (|TJ) . 
From these fits we obtain 1/e and the screening length due to free charges A^ defined as 
A^ 2 = A~ 2 — A^ 2 . These two quantities are the key quantities describing the KT charge 
unbinding transition; 1/e may be interpreted as describing the polarization due to bound 
dipole pairs whereas A^ can be interpreted as the Debye screening length related to the 
density rip of "free" charges A^ 2 = limpjeT . In the thermodynamic limit L oc A c — > oo 
all particles are bound into dipole pairs below the KT transition at T c whereas above T c 
some pairs are broken. This means that Xf = oo for T < T c and Xp < oo for T > T c . In 
accordance with this figure 2 shows how A^ 2 obtained in our simulations rapidly decreases 
as the KT transition is approached from above. Precisely at the KT transition one has the 
condition eT c = 1/4. |l]|| This is illustrated in figure 3 which shows e as a function of T 
for a sequence of constant particle densities. One notes that e increases monotonously with 
increasing T for low temperatures, goes through a maximum and then decreases towards 
e = 1 for higher T. Roughly this means that first the polarization due to bound pairs 
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increases because the average separation between the particles in a bound pair increases and 
then the polarization decreases because the number of bound pairs decreases due to thermal 
pair breaking at higher temperatures. The full curve in figure 3 corresponds to the condition 
IT = 1/4 and we use this as the determination of T c . This determination gives the phase 
transition line in the (n, T)-plane, as shown in the insert of figure 3. In the thermodynamic 
limit e has the critical behaviour 101 



e(T)-e(T c )oc±^/|T-T c | (19) 

where + and - refer to above and below T c . As seen in figure 3, the weak singular behaviour 
implied by Eq.(|19|) cannot be resolved by our present simulations. One notes, however, 
that the determined T c is close to the inflection point of the numerically obtained e-curve 
in accordance with Eq.(|19|). Associated with the weak singular behaviour of Eq.(plJ) is a 
corresponding singular behaviour of Xf || 

mA/oc-^=L= (20) 

as T c is approached from above. In figure 4 |lnA^ 2 | is plotted against 1/y/T — T c with 
T c determined from iT c = 1/4. As seen the critical behaviour given by Eq.(EUp is not 



discernible in the simulations. However, this result is expected because the true critical 
behaviour associated with Eq.(^) should in practice be extremely hard to resolve as a 
consequence of the extreme narrowness of the KT critical region. |17|,[18| In figure 4 we have 
also analyzed the data with respect to Eq. (^) following a commonly used procedure in the 
context of superconducting films and simulations on the 2D XY model: | In X F 2 \ is plotted 



against l/-y/T — T c where T c is a free parameter. As seen in figure 4 it is by this procedure 
possible to get a very good fit to Eq. (p0|) . Such fits are frequently claimed to be evidence for 
a KT critical behaviour. However, as discussed in ref. |L8| such fits do usually not reflect a 
critical KT property per se, but rather a property of the 2D Coulomb gas well outside the 
KT critical region. As is apparent from figure 4, our present simulations are consistent with 
this latter interpretation. 
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The low temperature phase displays a "quasi" 2D order in the sense that the correlations 
for large distances fall off like power laws. |L| In case of the charge density correlations we 
have 

sM = 0)cx-^y (21) 



for r » r where 



x{T) = ± (22) 



From a renormalization group (RG) point of view this means that each T < T c corresponds 
to a fixed point in the RG-flow. Q 

The RG-flow is towards vanishing density n so that for T < T c the line (n = 0, T) in the 
(n, T)-plane is a line of fixed points. Q Each such fixed point corresponds to a particular 
value of the critical index x(T). In figure 5 we show g(r, 0) as a function of r for a T below T c . 
The function g(r, 0) was obtain by directly measuring the charge correlations as a function 
of distance for the configurations generated by the simulation. The data for g(r, 0) is plotted 
as lng(r, 0) against lnr and according to Eq.([H]) the data should then fall on a straight line 
with slope x(T) for large r. As seen in figure 5 this prediction is borne out. The broken 
straight lines in figure 5 has the slopes given by 1/eT where e has been determined from 
l/e(k,u} = 0) as described in connection with figure 3. Thus the prediction x(T) = 1/eT 
is supported to high degree by our simulations. The fact that the power law decay of the 
correlations with distance and the power law index come out correctly gives us confidence 
in the present simulations. 

IV. /^-EXPONENT. 

Next we consider the non-linear response of the system when it is subject to an external 
force F ext = s^E where E is constant in space and time. This force generates a particle 
charge current I p . The charge current is in our simulations obtained from Eq.flHp. The 
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prediction is that below T c the generated charge particle current is a power law in the limit 
of small magnitudes of the force |J 

I P oc F: xt (23) 

In the context of a 2D superconductor the voltage V is proportional to the flux flow so 
that V oc I p whereas the force F ext is proportional to the Lorentz force so that F ext oc / 
where / is the external current applied to the superconductor. Thus in the context of a 
2D superconductor Eq.(p3|) corresponds to the non-linear JV^-characteristics for small / i.e. 
V oc I a . 

The question we are addressing with the present simulations is the value of the exponent 
a. There are two competing predictions: one is the AHNS-prediction |||6| 

a A HNS = 7^; + 1 (24) 
and the other is a scaling prediction fll2f 



a = — - 1 (25) 
IT y ' 

Figure 6 shows examples of the 7 p F ex . r characteristics obtained from our simulations. 
The data is plotted as ln/ p against \nF ext for a sequence of temperatures T. As is apparent 
from the figure the data fall to very good approximation on straight lines for small F ext , as 
predicted by Eq. (p3|) . The slopes of these lines give the values of the exponent a. The full 
curves in figure 6 are fits to the functional form 

I p = CF^e-^ 1 ^ 3 ^ (26) 

where a is the exponent and B and C are two constants. Fitting to this functional form 
turned out to be an expedient way of determining the exponent a: the a-values obtained by 
this fitting were the same as the ones obtained directly from the slope at small F ext but this 
latter procedure usually required much more computer time. 



11 



A heuristic motivation for Eq.(^) goes as follows: the particle current I p is proportional 
to the density of free particles np and the force F ext i.e. I p oc F ext np. The free particle den- 
sity may be related to a self-energy U se if for the creation of a free particle i.e. lnnp oc U se if. 
The self-energy corresponding to the effective interaction in Eq.(|TT|) is proportional to 
K (r /X) where the screening length A serves as an effective cut off of the particle interac- 
tion. |2| F~ xt has dimension of length and also serves as an effective cut off of the particle 
interaction. Consequently one may expect that whenever F~ x \ « A the effective cut off 
in the self-energy is proportional to F~ x \. This argument suggests that lnnp oc K (BF ext ) 
and Eq.(^) follows. 

Figure 7 shows the obtained values for the exponent a as a function of temperature T. 
These values are in the figure compared to the two competing predictions given Eqs (p4|) 
and (p5[), respectively. In this comparison we use the values of e obtained as described in 
section 2. As is demonstrated by figure 7, the scaling prediction given by Eq.(^5|) (full curve 
in the figure) is borne out to high precision whereas the AHNS prediction of Eq.(^) (broken 
curve in the figure) clearly disagrees with the data. One notes that the two predictions agree 
precisely at the temperature corresponding to 1/eT = 4 (crossing point between full and 
broken curve in figure 7). This corresponds to the critical condition for the KT-transition 
and to the universal jump value a = 3 at T c . |19|J5||| Above T c there are free charges even in 
the limit F ext = 0. Consequently one has I p oc n F (F ext = 0) ^ for very small F ext so that 
in principle a = 1 for T > T c . Thus in principle the exponents a jumps from 3 to 1 as T c is 
passed from below. However, as seen in figure 7, in practice the density of free charges np 
is dominated by the pair breaking mechanism also above T c for small F ext . This means that 
the exponent a corresponding to pair breaking above T c can in practice be determined to 
very good precision, as is apparent from from figure 7. From figure 7 we infer that the pair 
breaking exponent a is to very good approximation given by the scaling prediction Eq.(p5|) 
both below and above T c . 

The values of a given in figure 7 are for a fixed density n. In general the exponent a(T, n) 
is, of course, a function of both T and n. Thus we can also test the prediction for a as a 
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function of n. In figure 8 the data is plotted as a function of 1/eT for four different densities. 
The full straight line in figure 8 represents the scaling prediction of Eq. ( p5|) and the broken 
straight line the AHNS prediction of Eq.(^4|). As seen in figure 8 the data falls clearly on 
the full straight line representing the scaling prediction for all the various densities. Thus 
we conclude that the present simulations strongly supports the scaling prediction. 

An interesting observation in figure 8 is that the exponent a follows the scaling prediction 
(given by the full straight line) all the way down to a ~ 1 close to 1/eT « 2 at which point 
there is an abrupt crossover to a = 1. This suggests an abrupt crossover behaviour at 
T = 1/2 for small particle densities. The 2D Coulomb gas model has an equation of state 
which to leading order in the particle density n is given by p0[ 



p=(T-\)n for T> 1 - (27) 

and 

1 1 

p = -nT for T < - (28) 

where p is the pressure. For T < 1/2 this equation of state can be interpreted as the 
equation of state for an ideal gas of non-interacting dipole pairs. This suggests that the 
dominating part of the gas consists of dipole pairs in this small density limit. For such a gas 
of dipole pairs free charges can be generated by pair breaking caused by an external force. 
On the other hand for T > 1/2 the equation of state suggests a gas of free charges with no 
bound dipole pairs. This interpretation of the change of behaviour of the equation of state 
at T = 1 /2 is in accordance with the sharp crossover at T = 1/2 of the exponent a which is 
seen in figure 8. 

V. LARGE T-DEPENDENCE. 

In this section we focus on the large t-dependence of the charge density correlations 
below T c . Figure 9 shows our numerical data for the Fourier transform g(k,t). Our data 
suggest that the leading small k and large t-dependence is of the form 
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^2g— constk t 

g{k,t)<x (29) 

In order to establish this result we have in figure 9 plotted the logarithm of tg(k,t)/k 2 
against t for a sequence of fixed values of k. The form given by Eq.(|29|) implies that the 
data, when plotted in this way, should for large t fall on straight lines. Furthermore the 
slope of these lines should vanish as k approaches zero. As seen in figure 9 these features 
are very consistent with the data and the data for the smallest /c-values fall rather nearly 
on horizontal lines. The full straight lines in figure 9 are least square fits to the data in 
the region before too much noise sets in. For the four largest fc-values in figure 9 such 
lines can be determined without much uncertainty. The corresponding slopes, together with 
estimated uncertainties, are in figure 10 plotted against k 2 (filled circles with error bars). 
The broken straight line in figure 10 is a line through the origin which is least square fitted 
to the determined slopes. The fact that the slopes rather closely follow this line suggests 
that the slopes are proportional to k 2 for small k. The broken straight lines in figure 9 have 
the slopes given by the open circles in figure 10, i.e. they are the expected slopes for these 
smaller fc-values based on the /^-extrapolation of the slopes for the larger k- values. As seen 
in figure 9 the broken lines also fit rather well to the data, which lends further support to 
the conclusion that the slopes are indeed proportional to k 2 all the way down to k — 0. Thus 
we conclude that the data in figures 9 and 10 support that the small k and large t behaviour 
of g(k,t) to good approximation is given by Eq.(^). Simulations of the present type are 
of course always hampered by limited system sizes and time sequences. In particular we 
found that the smaller the fc-value the harder it was to obtain a large t-value free of finite 
size effects. E.g. the large t-part of the two smallest /c-values in figure 9 remain somewhat 
uncertain. Thus questions about logarithmic corrections to the leading t-dependence or 
non-leading terms appear to be outside the limitation of the present simulation precision. 



In ref. \12\ it was found that the function lim^o g(k, t)/k oc 1/t for large t in case of 



the 2D XY-model on a square lattice with TDGL (time-dependent Ginzburg-Landau type)- 



dynamics. This result was in ref. [|I^] associated with the vortex fluctuations. In the present 
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paper we confirm this conclusion by establishing the result directly in the Coulomb gas model 
with Langevin dynamics. In addition we have obtained the leading small fc-dependence for 
large t. 



comes from the small k. Thus we expect that the leading large t contribution to g(r,t) is 
given by 



for large t. The charge density is obviously a conserved quantity in the Coulomb gas model. 
Thus a pile up of charge in one place can only decay by diffusing away. Ordinary diffusion in 
2D leads to g(r, t) oc 1/t. However, from our simulations of the 2D Coulomb gas we conclude 
that the long range interaction between the particles changes this result into a more rapid 
decay g(r,t) oc 1/t 3 in the low temperature phase T < T c . In the high temperature phase 
T > T c the screening length A is always finite due to the presence of free charges. Thus in 
this case the decay of the charge density correlations are expected to decay exponentially, 
where the decay is dominated by a factor exp(— t\~ 2 const). We have not been able to 
explicitly verify this result in the present simulations, since the simulations are harder to 
converge in the high temperature phase. We note that, since in the small density limit the 
dipole pairs dominate the response for T < 1/2 (see the end of the preceding section), one 
might likewise expect that in practice the behaviour g(r,t) oc 1/t 3 , which we associate with 
the dipole pairs, also dominates the response in a region somewhat above T c (T c is always 
smaller than 1/4) for not too large time scales. Thus one might expect that the frequency 
response for small but not too small frequencies are dominated by the dipole pair response 
also in a region somewhat above T c . 

The scaling prediction for the exponent a given by Eq.(|25|) was in ref. [[12] based on 
the assumption that the charge density correlations g(r,t) can be associated with a scaling 
function A _2< I ) (rA _1 , t\~ z ) where z is the dynamical exponent and A is the screening length 
which diverges for any T below T c in the limit A c — > oo. Furthermore it was assumed that 



The form given by Eq. (|29|) implies that for large t the dominant contribution to g(r,t) 




(30) 



15 



the scaling function y) had the limits 0) oc x 2 ~ 1 / Te for large x and $(0,y) oc y^ 1 
for large y. Consequently we can now infer that the relation between g(r, t) and the invoked 
scaling form has to be 



\- z <5>(r\-\t\- z ) ocr 2 t 2 g{r,t) 



(31) 



since g(r ,t) oc t~ 3 for large t and g(r,t ) oc r 1 ' Te (compare discussion in connection with 
figure 5, to = Tq/D is the microscopic time scale and r is the microscopic size of a particle). 

From Eqs ( |12|) and ( |13|) one obtains the leading small a; dependence corresponding to 
Eq.([29l). Since according to Eq.(p9|) F(t) = limfc^o U(k)g(k, t) oc 1/t for large t, the leading 
small cu-dependences are proportional to 
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uj / dt sin Lut F(t) oc lu dx 
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respectively, 



pOO pl/oj 

uj dt cos ujt F(t) oc — uj I — oc uj In 
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and consequently 
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(32) 



(33) 



(34) 



Im 



const uj In a; 



(35) 



for small uj. The constant const depends on T and has the same value in the last two 
equations. 

It is instructive to translate this result into the linear response function corresponding 
to the conductivity o~{uj) for a 2D superconductor. The connection is cr{uj) oc [— iue(0, 
Consequently we have the leading small (^-dependence given by 



i-i 



Rea{uj) oc — In \u\ 



(36) 



and 
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or equivalently 



constant n 
Ima{ui) oc h — sign(o;) 



constant . . s 

o\uj) oc ln(— iuj) 

—iuj 



(37) 



(38) 



For a 2D superconductor the dissipation for small u are due to thermally created vortices. 
Hence Eq.([36|) suggests that for a 2D superconductor the real part of the conductivity 
diverges logarithmically below T c . 

VI. PHENOMENOLOGICAL DESCRIPTION. 

As described above, we compute the Fourier transform of the charge density correla- 



tion function g(k,t) in our simulations. By aid of Eqs ( |12l) and ( |13l) we obtain the fre- 
quency dependence of the dielectric function l/e(k,u). In figure 11 we present the result for 
Re [l/e(k, uj)] and Im [l/e(k, uj)] for the smallest k we managed to converge (k = 0.039 (r^ 1 )). 
As is apparent from figure 9 we expect that for such a small value of k the real part of 1/e 
should vanish linearly with u, as discussed in connection with Eq . (p4|) . For large uj, on the 
other hand, the dipole pairs have no time to respond so that in this limit we expect that 
Re [l/e(k,uj)] = 1. Thus we expect the real part for small enough k to be of the form 
Re [l/e(Jfe,u;)] « Re [l/e(0, uj)} where 



Re 



e(0,uj) 



e(0,0) 



+ 



1 



e(0,0) 



UJ 



uj + G(uj) 



(39) 



provided < G(u) < oo and G(u = 0) = constant > 0. If G(u) only depends weakly on 
uj, we can approximate G(uj) by a positive constant G(uj) ~ Uq. Using this approximation 
we can now obtain the corresponding approximation for the imaginary part by using the 
Kramers-Kronig relation. This leads to 



Re 



1 



1 

e(0,0) 



+ 



1 

e(0,0)_ 



UJ 



UJ + LU 



(40) 



and 
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One notes that Eqs fl40|) and (41) correctly reduces to Eqs (|34]) and (|35| ) in the small uj- 



limit. The form of the a;-dependence given by Eqs (|40|) and (|T) is identical to the form given 
by the MP-description. The MP-description was originally motivated from a heuristical 
argument for the dipole pair response. As described here, we can also view the MP-form 
as a simple interpolation between two known limits. In figure 11 we have fitted Eqs ([40"D and 
([|l]) to the data from the simulations by using the two constants ujq and 1 — as fitting 

parameters. As seen from the figure, the data is very well described by the MP-form. From 
the fitting we obtain l/e(0, 0) = 0.91. This is fairly close to the value obtained directly from 
the simulations l/e(k, 0) = 1/e ~ 0.92 (compare figure 9, k~ 2 g(k, t = 0) oc l — l/e(k, uj = 0)). 
The value obtained for ujq from the fitting was ujq ~ 0.36 ± 0.02. This can be compared to 
G(0) ~ 0.24 which is obtained directly from the simulations: From the data in figure 9 we 
get g(k,t)/k 2 = A/t for small k and large t. By aid of Eq.(p2|) we then obtain the value of 
the constant in front of the linear small u dependence. Finally, by using l/e(k, 0) from figure 
1, we obtain G(0) as G(0) ~ (1 — l/e(k,0))/ constant. One notes that G(0) is somewhat 
smaller but of the same magnitude as the ujq determined from the fitting to the data, which 
supports the assumption that G{ui) only has a quite weak tu-dependence. Thus we conclude 
that Eqs (|40] and flHf) , obtained from Eq.(^) by using the approximation G(oj) Pd ujq, gives 
a very good and quite consistent description of the data from our present simulations. 

So far we have focused on the cu-dependence in the small fc-limit. Next we consider how 
this u dependence is changed as k is increased. Figure 12 shows the real and imaginary part 
of l/e(k,u) as a function of to for a sequence of increasing k- values. As seen in figure 12 
the imaginary part appears to be remarkably independent of k over the range of /c-vales in 
the figure (0.039 < k < 0.31), whereas the real part decreases significantly with increasing 
k. However, since the real and imaginary part are related by a Kramer-Kronig relation, 
the change in the real part part must have a corresponding change in the imaginary part. 
The point is that the corresponding change in the imaginary part is concentrated to small 



frequencies. Since the imaginary part is almost independent of k whereas the real part 
changes significantly, one concludes that Eqs ([It]) and ( |41| ) does not describe the data as 
k is increased, as expected from the motivation of these equations. In order to quantify 
this feature in a practical way we focus on the change with k precisely at the maximum 
of the imaginary part. As seen in figure 12 this frequency is close to u ~ 0.36 and is 
approximately independent of k over the k range in the figure . Now the MP-form given by 
Eqs fl40D and ([41~1) has the property that, at the maximum frequency for the absolute value 
of the imaginary part the, ratio between the absolute value of the imaginary and real part 
is precisely 2/7T. In figure 13 we have plotted this ratio as a function of k. For small k this 
ratio approaches 2/ir as expected from Eqs ( [40|) and (|4TD. However as k increases this ratio 
also increases. Ordinary diffusion corresponds to 



1 1 

+ 



e(0,u) e(0,0) 



1 

1 - 



(42) 



uo + iDk 2 



e(0,0) 

The real and imaginary part is for this case of the usual Drude form. The ratio between the 
absolute value of the imaginary and real part precisely at the maximum of the absolute value 
of the imaginary part is in the Drude case unity. It appears from figure 13 that the ratio 
approaches the ordinary Drude value unity as k increases. This change of behaviour from MP 
to Drude is somewhat reminiscent of Eq. (|29|) where the diffusion like factor exp(— const k 2 t) 
becomes more important with increasing k. 

VII. CONCLUDING REMARKS. 

The model investigated is a model of 2D vortex fluctuations. Consequently, the prop- 
erties of this model should be directly reflected in measurements on 2D superconductors 
like superconducting films and 2D Josephson junction arrays. In our simulations we verified 
that the KT critical region is very narrow. For the resistance R of a superconductor the 
predicted KT-critical behaviour is || 

In R oc 1 (43) 



' c 
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where T CG is the effective temperature variable for the vortices. The narrowness of the 
critical region means that in practice the KT-critical behaviour probably cannot be resolved 
in resistance measurements, as has been pointed out earlier. P, P^7| , P^] In our simulations we 
obtained T c to good precision. However, if we instead used T c as a free parameter, then we 
showed that the data could indeed be fitted to the KT-critical behaviour over several decades. 
This procedure of using T c cis ci free parameter and fitting to Eq. (fE|) is the commonly used 
way of establishing the KT-critical behaviour for the resistance R of 2D superconductors. It 
is frequently found that by treating T c as a free parameter the data can be nicely fitted to 
Eq . (f43"D , as was also the case for our simulations. However, in our simulations we also found 
that the T c from the fit was significantly different from ( i.e. 13% lower than) the true T c . 
Thus we conclude that these type of fits have in fact no direct bearing on the "real" critical 



behaviour. |fr?| , p!8| It should also be noted that although the fitted T c for a superconductor 
can appear to be quite close to the true T c in real temperatures, the corresponding difference 
in the effective vortex temperature variable T CG is usually much larger. The crucial point 
here is that our simulations suggest that, if the true T GG is used, then the data does not 
follow the functional form given by Eq. (pf). 

We also obtained the non-linear TV-exponent a from the simulations and verified that 



1/eT — 1 plf , as proposed in ref. |l2j and which is quite different from the earlier 



AHNS-prediction |5||| clahns- However, precisely at T c both prediction give a = 3 since 
a = 2(1ahns ~ 3- Below T c the new value is larger than the AHNS-value. In principle a = 1 
above T c due to usual flux flow resistance of free vortices. Nonetheless, we found that in 
practice a non-linear TV-exponent describing the pair breaking could also be determined to 
very good precision above T c up to roughly 2T C at which point there was a rapid cross over 
to a = 1. All the way up to roughly 2T C we found a = 1/eT — 1 to very good approximation. 

For a 2D superconductor the exponent a is directly related to l/eT GG oc p(T)/T where 
the proportionality constant is a combination of fundamental physical constants and p(T) is 
the 2D superfluid density. Thus in order to test the prediction for a one needs to know the 
temperature dependence of p(T). One way is to measure the complex impedance Z(T,u) 
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since 



rpCG 

Z(T, oj) = -iwL k {T)e(k = 0, u, T) oc -^-e(fc = 0, u, T) (44) 

where Lk(T) is the sheet kinetic inductance and the proportionality factor is again just a 
combination of fundamental physical constants. Consequently, if Z(T,u) is measured for 
very small frequencies for a sample and the non-linear IV-characteristics for the same sample 
is measured to high precision, then the prediction for a for T < T c can be directly put to 
experimental test. Alternatively, if Z(T,u) is measured for a small but finite frequencies, 



then the u = 0-limit can be extracted by using Eqs (f40|) and (fy]). For T > T c one needs 
both Lk(T) and e(T). The sheet kinetic inductance Lk(T) can often be quite well determined 
from the complex impedance. The qualitative behaviour of e(T) is clear from figure 3. It is 
also possible to make a somewhat more quantitative determination of e above T c by starting 



from e(k = 0,u,T) at a finite small frequency and fitting to the MP- phenomenology. p2|JTT 
We concluded from our model that 2D vortex fluctuations has a long range 1/t-tail in the 
vortex correlations below T c . This means that the conductivity <j{u) is of the form given by 
Eq. (|55D . The crucial feature is that the next leading term for small u is a logarithm ln(— iu). 
Such a logarithmic term is strongly supported by experiments on a 2D array of Josephson 
junctions. |§ These experiments are in fact very well described by the MP-phenomenology. 
@ In connection with these experiments there has been other proposal for the origin of this 



logarithm, e.g like vortex-spin wave coupling [£3] or as a specific single vortex property of 



a proximity coupled array 24]. The point we are making here is that this logarithm is an 
intrinsic collective property of 2D vortex fluctuations which is strongly linked to the long 
range logarithmic vortex interaction. 

The general form of the frequency dependence is given by Eq.(p9|) which reduces to the 
MP-form of Eqs (|40| ) and fl4I| ) for G(ou) = const. The point here is that, since there is 
no a priori characteristic frequency scale other than the microscopic t$ 1 one expects that 
G(u>) ~ G(0) as long as u « . We believe that the fact that the MP-form describes a 



variety of experimental data as well as simulations very well [H,|?Hll|li reflects this aspect: the 
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MP-form describes the dynamics for frequencies much smaller than the basic microscopic 
frequency scale and the functional form is independent of the microscopic details of the 
dynamics. In this sense it describes a universal behaviour of 2D vortex fluctuations. 

One may also note that Uq(T) in Eqs ( f4"0"|) and ( |4TD depend on T. Thus the MP-form 
can be tested either for fixed T by varying u or for fixed u by varying T. This latter way is 
more common in experiments. A particular feature of the MP-form is the fact that the ratio 
between the imaginary and real part of the response is 2/ir precisely at the dissipation peak. 
The position of this peak is often quite well defined in the experiments and peak ratios close 
to 2/tt have been found both in experiments and simulations. [f7HTTH In the present paper 



we directly verify that the peak ratio 2tt is an intrinsic dynamical small frequency property 
of 2D vortex fluctuations. 
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FIGURES 

FIG. 1. The static dielectric response function l/e(k) plotted as a function of the wave vector k 
(in units of2ir/L where L is the system size ) at a small constant particle density, n = 0.005 (tq 2 ). 
The filled circles represent the data from simulations at T = 0.12, 0.14, 0.16, 0.18, 0.20, 0.23, 
0.26, 0.29, 0.32, 0.35, 0.38, 0.41, 0.50, from top to bottom, respectively The full curves are fits of 
equation 

FIG. 2. The inverse square of the screening length X F 2 extracted from the data in figure 1 by 



using Eq.(18). The open circles represent Xp 2 (in units of (2tt/L) 2 and n = 0.005 (r^ 2 )) plotted 
as a function of temperature T. The critical temperature T c is also shown ( T c is determined 
as described in connection with figure 3). The figure illustrates the charge unbinding transition: 
Below T c the charges are bound together in dipole pairs. These pairs start to unbind at T c , resulting 
in a rapid increase in the density of free charges tx A^ 2 as T passes T c from below. 



FIG. 3. The dielectric constant e extracted from the data in figure 1 by using Eq.(lt). The 



various symbols represent, from bottom to top, the determined value of l/e(k) for the particle 
densities, n = 0.001, 0.005, 0.01, 0.025, 0.05, 0.075, 0.025, 0.05, 0.1 (r^ 2 ), respectively. The broken 
curves are guides to the eye. The full curve represent the KT critical condition e = 1/4T C . The 
full curve crosses the broken curves close to the inflection points of the broken curves as expected 
for the KT-transition. The crossing points between the full curve and the broken curves gives the 
phase transition line in the (n, T)-plane as shown in the insert. The estimated error bars in the 
insert are due to numerical uncertainties as well as finite size effects. It should be noted that the 
uncertainties in the values of T c are very small. 
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FIG. 4. Test of KT-criticality: The quantity | In A^ 2 | is plotted against l/y/T-T c . Plotted 



in this way the data should, according to Eq.(2C), fall on a straight line, provided the data is in 



the critical region. The data for Ap is the same as in Egure 2. The filled circles represent the case 
when the true T c is used (T c = 0.215(7), determined as described in connection with hgure 3). The 
broken curve is a guide to the eye. As expected Eq. $2(\) does not describe the data because of the 
unusual narrowness of the KT critical region. The critical behaviour would only show up closer to 
T c for much larger values of Xf which cannot be converged in the present simulations. If instead 
T c is treated as a free parameter then data can indeed be manipulated to fall on a straight line. 
The open circles show the case for T c = 0.187. In this case the data fall on a straight line over 
several decades as indicated by the dotted straight line in the figure. 

FIG. 5. The charge density correlation function g(r, t = 0) as a function of distance r for 
T = 0.18, n = 0.005 (circles), and T = 0.16, n = 0.025 (diamonds), respectively. The data is 
plotted as |logg(r, 0)| against logr in order to test the prediction g(r,0) oc r _1 / eT for large r. 
Plotted in this way the data should fall on straight lines for large r and this prediction is borne 
out. The broken lines have the slopes given by 1/eT where e has been determined as discussed in 
connection with Egure 3 and fit the data very well. Consequently the prediction g(r, 0) oc r~ l / eT 
is verified to high degree by the present simulations. 

FIG. 6. The charge current I p as a function of an external force F ex t = S{E for the temperatures 
T = 0.12, 0.14, 0.16, 0.18, 0.23, 0.26, 0.29, 0.50 at constant density n = 0.005 (r^ 2 ). The data 
from the simulation are given by the filled circles with error bars and is plotted as log I p against 
logF ext . The full curves are fits to Eq.p6|). From these fits the exponents a dehned by I p oc F£ xt 
are determined. Alternatively the exponent a can be determined directly from the slopes at small 
F ext (see text). 
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FIG. 7. The non-linear IV-exponent a obtained from the data in figure 6. The exponent a is 
plotted as a function of temperature (filled circles with error bars). The full curve is the scaling 



prediction Eq.(25) and the broken curve is the AHNS prediction Eq.($4). The values ofe needed to 
make the comparison were determined as described in connection with figure 3. The data strongly 
favors the scaling prediction. 

FIG. 8. The IV-exponent obtained from the simulations at four different particle densities 
(n = 0.001, 0.005, 0.01, and 0.025 (r^ 2 )). The data is plotted as a against 1/eT. Plotted in this 
way the scaling prediction corresponds to the full straight line and the AHNS prediction to the 
broken straight line. The vertical error bars estimate the uncertainty in the value a for a given T 
and the horizontal the uncertainty in e. The data clearly verifies the scaling prediction. Precisely 
at T c one has 1/eT = 4 and both predictions give a = 3. Note that a in practice can be determined 
also above T c (see text) and clearly follows the scaling prediction all the way upto 1 /eT ~ 2 at 
which point there an abrupt cross over to a = 1. 

FIG. 9. The time dependence of the charge density correlation function g below T c (T = 0.18 
and n = 0.005 (tq 2 )). In order to verify the t-dependence g(k,t) oc k 2 e~ constk2t /t given by 
Eq.(2£) the data is plotted as log(tg(k,t)/k 2 ) against t. The data for large t should then 
fall on straight lines where the slopes go towards zero as k is decreased. The data is given 
by the open circles and the eight data sets correspond from top to bottom to the k-values 
k = (1,1), (2,2), (3,3), (4, 4), (5,5), (6,6), (7, 7) and (8,8) (in units of 2ir/L). The t-dependence 
of Eq.(^) is apparently borne out to fairly good approximation. The full straight lines are fits to 
the data for the larger k-values before too much noise sets in. The broken straight lines have slopes 
extrapolated from the full straight lines using Eq.ffijQ as described in connection with figure 10. 
Also the broken straight lines fit the data fairly well. 
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FIG. 10. The k-dependence of the slopes in figure 9. The slopes of the full straight lines in 
figure 9 are plotted against k 2 (filled circles with error bars). According to Eq.^9^) these slopes 
should extrapolate linearly to zero. The broken straight line is a line through zero which is fitted 
to the filled circles and shows that the slopes apparently to good approximation are proportional 
to k 2 . The open circles give the expected slopes for some smaller k-values and correspond to the 
broken straight lines in figure 9. 

FIG. 11. The real (circles) and imaginary (diamonds) part of the frequency dependent di- 
electric constant l/i(k,u>) for a small k (the h-vector is (1, 1) (in units of2n/L) and the data is 
the same as the top data set in figure 9). The full drawn curves are fits to the MP-form given by 



Eqs (|4Q ) and (41) using the two constants ujq and 1 — l/e(0, 0) as free parameters. The MP-form 



describes the data very well. 

FIG. 12. Real and imaginary parts of l/e(k,uj) for the same parameters as in hgure 9. The 
real part of l/e(k,uj) is represented by the upper set of curves. The amplitude decreases with 
increasing k (the k-values are the same as in figure 9). The uppermost curve is the same data as 
was shown to be well described by the MP-form in figure 11. The k-dependence of the imaginary 
part (lower set of curves) is by contrast very small. Note that the absolute value of the imaginary 
part has a maximum at uj ~ 0.36 for all k-values in the figure. The hgure illustrates that the 
MP-form as expected only describes the data in the limit of small k. 

FIG. 13. The peak ratio of the imaginary and real parts of l/e(k,u) as a function of k. The 
peak ratios are obtained from the data shown in hgure 12 and are denoted by open circles. The 
broken curve is a guide to the eye. For small k the ratio approaches the predicted MP value of2/n 
whereas the Drude ratio of unity is approached for large k. 
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